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Abstract 

We present a Monte Carlo study of the bond and the site directed percolation model on the 
(d+ l)-dimensional simple-cubic and the body-centered-cubic lattices with d = 2 to 7, as well as on 
the square, triangular, honeycomb, and kagome lattices for d = 1. With a cumulative-probability 
and a linking-hashing technique, the simulation is reasonably efficient, independent of d, and can 
be up to very large virtual lattices (e.g., more than 10 27 sites for d = 7). A dimensionless ratio 
Qj\f is defined to locate the percolation threshold p c . Universal scaling of Qjy near p c is checked, 
and accordingly the Qj\f data for systems in the same dimension are simultaneously analyzed by 
a formula, in which universal parameters occur only once. In comparison with existing results, 
our estimates for p c are comparable or somewhat better for d = 1, and much more precise for 
d > 1. Extensive simulations are then carried out at the estimated thresholds. The probability 
distribution and the finite-size scaling of several quantities are studied. The critical exponents are 
accurately determined. At the upper critical dimensionality d c = 4, logarithmic corrections are 
clearly observed and found to be well described by the theoretical predictions. For d > d c , the 
mean-field behavior is confirmed. 
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I. INTRODUCTION 



Directed percolation (DP) is a geometric model generalized from ordinary (isotropic) 
percolation, first introduced by Broadbent and Hammersley (1957) 1). It is a generalized 
percolation with a special direction along which activity can only propagate in one way but 
not the reverse lrM|. This special direction is normally considered as the temporal direction 
while the others are spatial. It is widely agreed that a variety of phenomena in nature are 
related to DP, including epidemic disease |2|, forest fire [l[ O], transport in porous media, 
etc. |3, 4|. 

As a typical non-equilibrium process, DP has attracted much attention from statistical 
physicists because of its simplicity and rich behavior 0Q. Unfortunately, analytical so- 
lutions are scarce for DP, even in (1 + 1) dimensions. The direct experimental explorations 



are also difficult; until now^ 



turbulent liquid crystals 



;he exclusive example to our knowledge is the DP criticality in 
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[x3 1 . This makes numerical means indepensible to study the 
DP process, which include series expansion, transfer matrix, Monte Carlo simulation, etc. 
Fruitful results have been obtained. By series expansion, Blease provided the first estimate 
of the percolation threshold p c for bond DP on most of the typical lattices in (d + 1) di- 
mensions with d = 1-7 jl4|. Essam and Jensen et al. gradually improved the precision for 



some (l+l)-dimensional lattices, both for bond and site DPs [15M22I]. Similar progresses 
have been achieved in (2+1) and higher dimensions. Grassberger and Adler et al. improved 
the precision or established the first estimate of p c for several lattices 23 2s| . Among the 
existing results, the estimate of p c on most of the (l+l)-dimensional lattices is impressively 
accurate; the uncertainty is at the 8th decimal place. An exception is the site DP on the 
honeycomb lattice 18j . of which the error margin is at the Qth decimal place. In (2+1) and 



higher dimensions, the precision of p c is less satisfactory. 

Numerical methods like Monte Carlo simulation are usually limited to finite systems. 
Thus, theory of finite-size scaling (FSS) is widely used to extrapolate to infinite systems. 
Underlying FSS theory is the hypothesis of universality. It states that the singular critical 
behavior for various seemingly unrelated systems is governed by the same set of critical 
exponents. In equilibrium statistical mechanics, the universality hypothesis is supported 
by renormalization group (RG) theory. In this framework, critical systems of the same 
universality are in the same critical surface in the parameter space constructed by scaling 



fields, and their scaling behavior is governed by a common fixed point. In two dimen- 
sions, beside numerical results, evidence for the universality hypothesis further arises from 
analytical calculations, conformal field theory, and Coulomb-gas theory Compared to equi- 
librium statistical mechanics, less is known for non-equilibrium phase transition. This is 
due to the absence of an analog to equilibrium free energy and thus of a general treatment. 
Nevertheless, Janssen and Grassberger did propose a well-known universality hypothesis for 



DP 
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11] : short-range interacting models, which exhibit a continuous phase transition into 
a unique absorbing state, belong to the DP universality class, if they have a one-component 
order parameter and no additional symmetries. 

RG theory not only suggests the same set of critical exponents for systems in a university 
class, but also implies universal scaling functions. In FSS theory, the linear system size is 
assumed to act as a scaling field with renormalization exponent —1, and thus FSS functions 
are also universal. Universal FSS behavior has been carefully checked for standard perco- 
lation and the Ising model, two paradigmatic systems in equilibrium statistical mechanics. 
However, for DP, less is known. 

The FSS behavior of observables, such as specific heat and susceptibility in equilibrium 
system, is typically described by divergence according to power law modified by corrections 
to scaling. Thus, numerical treatment of the Monte Carlo data has to include both unknown 
critical exponents and unknown corrections. This introduces complications for data analysis 
and also restricts the accuracy of results. To overcome this difficulty, Binder introduced a 



dimensionless ratio = 1 — (m 4 ) /3(m 2 ) 2 for t 
density and ( ) represents statistical average 



revising model, where m is the magnetization 
From universal FSS function, it can be 



shown that, if corrections to scaling are ignored, the U4 values for different system sizes are 
identical to a common value at the critical point, and further this common value is universal. 
The U4 data for different sizes typically display an approximately common intersection, 
which can be used to locate the critical point. Analogous universal ratios are widely used in 
systems besides the Ising model, and have been proved to be very powerful in numerically 
determining transition point js0|. 

In this work, we shall define a dimensionless ratio Qa/ - for DP on the basis of the number 
of wet sites Af '. It is observed that, as expected, Qj\f does obey universal FSS behavior for 
systems in the same universality. Further, as for equilibrium system, the dimensionless ratio 



Qjsf proves to be very useful in determining p c . We then follow the strategy in Ref. 3l| and 



carry out a simultaneous analysis of systems in the same universality class, which include 
all site and bond DPs in the same dimension d. In such an analysis, all the Monte Carlo 
data with the same dimensionality d are combined together, and universal parameters occur 



only once in the fitting formula. This further improves our estimate of p c . 



reason of the simultaneous analysis has been discussed in detail in Ref. 31], arid shall be 



'he underlying 



briefly mentioned. Finite-size data usually suffer from corrections to scaling. In the first 
order, these corrections are proportional to the leading irrelevant renormalization field, of 
which the amplitude approximately parameterizes the distance between the fixed point and 
the location of the model in the critical surface. Different systems in the same universality 
class usually have different amplitudes. A model with small irrelevant field can be well used 
to determine universal parameters like the critical ratio, but is not suitable to determine the 
irrelevant exponent. The well estimated universal ratio can help the Monte Carlo data for a 
model with a relatively large correction to determine the irrelevant exponent. Accordingly, 
the estimate of percolation threshold and other universal parameters is improved. 

Extensive simulations were then carried out at the estimated percolation thresholds, 
and universal FSS functions of several quantities are studied. These include the wetting 
probability V, the number of wet sites A/", and the mean gyration radius 1Z. To study their 
FSS functions, we measure the probability distribution (histogram) of these quantities. The 
critical exponents are determined, and the (hyper-) scaling relations are checked for d < d c , 
where d c = 4 is the upper critical dimensionality. Logarithmic corrections at d — d c are 
clearly observed and analyzed according to the theoretical calculation. For d > d c , the 
hyperscaling relation breaks down and the mean-field prediction is confirmed. 

The remainder of this paper is organized as follows. Section ILT1 briefly introduces the DP 
model, defines the lattices, and describes the algorithm. Results are presented in Sees. II III 
and IIVt which include the determination of the percolation thresholds and the study of the 
critical FSS behavior. A brief discussion is given in Sec. |Vj 
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(a)square (b)triangular (c)honcycomb (d)kagome 

FIG. 1: Some (l+l)-dimensional lattices. The horizontal and vertical directions are regarded as 
the spatial and temporal directions, respectively. 



II. MODEL, LATTICE, AND ALGORITHM 



A. Model 



Among various definitions of DP, the model of growing from a single wet site and the rule 
of parallel updating are adopted in this work. Consider the square lattice in Fig. [Ha), where 
the horizontal (vertical) direction is the spatial (temporal) direction. One starts with a single 
wet site at time t — 0. At each time step t, in the bond DP process, every edge, connecting 
an already- wet site at t and a lattice site at t + 1, becomes occupied with probability p; sites 
at t + 1 with one or more occupied edges are said to be wetted. In the site DP process, at 
each time step t, every lattice site at time t+1, with one or more wet neighboring sites at t, 
is turned wetted with probability p (sites with no wet neighboring site at t cannot become 
wetted). For occupation probability p > 0, all the lattice sites in Fig. HJa) can potentially 
be wetted, while those outside the geometry have no chance to get wetted. 

For d > 1, there exist two distinct phases in the DP process. For small p, the number of 
wet sites J\f(t) rapidly approaches to as time t elapses, i.e., the propagation quickly dies 
out. For large p, the DP process becomes absorbed and J\f(t) increases linearly as t d . A 
phase transition point p c occurs to separate these two phases. The number of wet sites Af(t) 
behaves as 



e -*/£ii(p) p<p, 



c 



Af(t) oc < 



t e p = p c (!) 

t d [a + 0(e-* /€ i: (p) )] p > p c 

where < a < 1 accounts for the fraction of wet sites over the system. The correlation 
length £|| in the temporal direction is a function of p, and diverges as oc \p — p c \~ Un when 
p — > p c . The critical exponent < 9 < 1 characterizes the scaling behavior of N"(t) at p c . 



(a)SC (b)BCC 
FIG. 2: (2+l)-dimensional SC and BCC lattices. 



B. Lattice 

In (1 + 1) dimensions, besides the square lattice, we also study bond and site DPs on the 
triangular, honeycomb, and kagome lattices. On the triangular lattice (Fig. QJb)), at each 
time step t, the propagation of wet sites can be from time ttot + 1 and from t to t + 2. On 
the honeycomb (Fig. QJc)) and the kagome lattice (Fig. [T^d)), the unit of temporal direction 
is defined such that the propagation occurs between time t and t + 1/2 and between t + 1/2 
and t + 1 . 

For d = 2 to 7, we study DP on the simple-cubic (SC) and the body-centered-cubic 
(BCC) lattice. On the SC lattice, each lattice site at time t has A = d + 1 neighboring sites 
at t + 1; on the BCC lattice, one has A = 2 d . Figure |2] shows the SC and BCC lattices in 
(2 + 1) dimensions. The visualization becomes difficult for d > 2, and we shall define the 
SC and the BCC lattice according to the association between a site at t and its neighboring 
sites at t + 1. This will be used in the actual Monte Carlo simulation. In (d+ 1) dimensions, 
let v t := (xi, x 2 , ■ ■ ■ , Xd) specify a lattice site at time t, where integer xi > (i — 1, 2, . . . , d) 
is the spatial coordinate. On the SC lattice, the coordinates of the A = d + 1 neighboring 

sites at t + 1 are (xi, £2, • • • , x<i), (xi + 1, £2, • • • , x<i), , (#1, £2, • • • , %d + 1), respectively. 

In other words, one of the A = d + 1 neighboring sites at time t + 1 of v< has the same 
spatial coordinates as v t , while the others differ from v t by one lattice unit in one of the d 
directions. On the BCC lattice, the spatial coordinates of the A = 2 d neighboring sites at 
time t + 1 of Vf read (x\ ± 1, x 2 ± 1, . . . , ± 1). Namely, the coordinate in each of the d 
spatial directions either increases or decreases by one lattice unit. 



C. Algorithm 



As mentioned above, all lattice sites within the geometry in Figs. [T] and [2] have non-zero 
probability to get wetted. Thus, the number of the potentially wetted lattice sites V(t) 
increases approximately as a(t + as t elapses, with a > a constant. On the square, 
triangular, honeycomb, and kagome lattices, constant a is 1/2, 1/2, 1, and 3/2, respectively; 
for the (d + l)-dimensional SC and BCC lattices, one has a ~ l/(d + 1)! and l/{d + 1), 
respectively. Recording all these lattice sites would request a huge amount of computer 
memory for large t and d. Fortunately, this can be avoided in the DP process. When 
propagating from ttot + 1, one only needs to know the locations of the wetted sites at t; the 
unwetted sites are irrelevant in the propagation procedure. Thus, we introduce two first-in- 
first-out queues Q t and Qt+i to store the coordinates of the already- wet sites at t and of the 
unfinished list of wetted sites at t + 1 (for the triangular lattice, an additional queue Qt+2 
is needed), respectively. Each element of Q is a d-dimensional vector v := (xi,X2, ■ ■ ■ ,Xd), 
specifying the spatial coordinates of a lattice site. With these queues, the lattices in Figs. [U 
and P2] are not directly stored in computer memory, and thus become "virtual" . 

For large d, the percolation threshold p c approximately equals to 1/A, and thus p c is small 
for large d. Therefore, near the critical point, most of the lattice sites are unwetted, and 
storing the wetted sites instead of the whole virtual lattice would save numerous computer 
memory. 

The small p c also implies that many random numbers and many operations would be 
needed if all the potentially occupied edges (sites) are visited and determined whether or 
not they are occupied. The simulation efficiency would drop quickly as the coordination 



number A increases. A more efficient procedure follows 
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Let Vt be a wet site at t, 
one counts its A neighboring sites sequentially, and defines P{k) := (1 — p)^" 1 ^ to be the 
probability that the first (k — 1) neighboring sites are empty (unwetted) while the kth site 
is wetted. The cumulative probability distribution C{k) is then 

k 

C(k) = Y,P(k') = l-(l- P ) k ■ (2) 

k'=l 

Suppose the current wetted site is the fcoth neighbor of v t , one can obtain the next to-be- 
wetted site as the (ko + k)th neighbor of v t by drawing a uniformly distributed random 



number < r < 1 and requesting 

C(k - 1) < r < C(k) . (3) 

The value of k can be obtained within log 2 (A) steps by the standard bisection search method, 
which compares r to the value of the middle element of C(k) in the searching interval. 
Thus, instead of visiting all the potentially wetted sites, one can directly jump to the next 
to-be- wetted site. This step is repeated until all the neighboring sites of v t are visited-i.e., 
ko + k > A. By this way, the simulation efficiency is significantly improved for large d. 
On this basis, one can formulate the algorithm as 

A. Visit sites in queue Q f sequentially. For each site w t in Q f , starting with the fcoth neighbor 

of v t (ko = at the beginning), 

A.l Calculate the location (fc + k) of the next to-be-wetted site v t+1 by Eq. 03]). 

A. 2 Check if v t+ i has already been wetted (visited) for bond (site) DP. If not, make 
v t+1 wetted, and update queue Qt+i- 

A. 3 Set ko -f- ko + k. Repeat A.l and A. 2 until all the A neighboring sites of v t are 
visited. 

B. Repeat step A until all the elements in Q t are visited. Set Q t <— Qt+i- 

The simulation continues until the aimed maximum time step t max is reached or the prop- 
agation dies out. In our Monte Carlo simulation, the maximum time t max is set in range 
2 13 < t < 2 16 

In the actual implementation, some difficulties may occur in step A. 2. In this step, for 
bond DP, the main task is to check whether a to-be- wetted site (derived from Eq. ([3])) is 
already wetted. For d=l and 2, this can be done easily by introducing an auxiliary "flag" 
to specify the wetted sites, since the memory is not yet a problem in modern computers 
typically of 2 GBytes or more. For d > 2, especially for large d, storing all lattice sites 
becomes impractical. In this cs.sc, clS mentioned earlier, one introduces a to-be-complete 
queue Qt+i to store the wet sites only. When judging whether or not a to-be- wetted site 
v t +i should be added to the queue Qt+i, one may compare the coordinates of each already- 
in-Q t+1 site to those of v t+1 , one by one. Such a comparison need about (N(t) 2 ) operations, 

o 



FIG. 3: Sketch of the linking-hashing technique. All elements in the same linked class have the 
same first two coordinates (x\,X2). The linking direction is backward in queue Q. 

with N(t) is the number of wetted sites at time t. This effectively leads to a computational 
critical slowing-down. 

We use a linking-hashing trick to resolve this difficulty. One classifies the sites in Qt+i by 
their first two coordinates (xi,X2), and links all the sites with the same (xx,^) by directed 
pointers. This is shown in Fig [3j where the pointer direction is backward (the sites at 
the head of the arrow are added to queue Qt+i earlier than the sites at tail). A hashing 
table M(xx,X2) is introduced to store the top element of each class. With these auxiliary 
informations, the status-checking procedure can be done as following. When v 4+1 is derived, 
the linked class is immediately identified by its first two coordinates (aq,^), and the top 
element of the class is read out from table H(aq,X2)- Then, along the directed linking, it is 
checked whether Vf+i is identical to any site in this class. For completeness, we re-express 
step A. 2 in the aforementioned algorithm as 'A. 2 Check if v t+1 was already wetted (visited) 
for bond (site) DP. If not, make v t+1 wetted, and update queue Qt+i, the linking, and the 
hashing table'. 

At the percolation threshold p c , the average number of wet sites scales as J\f{t) oc t e . The 
scaling exponent 9 is in range 0(d > 4) < 9 < 0.313686(d = I). More precisely, one has 
N(t) ~ I for d > 4-i.e., only a few wetted sites survive on the average. This means that 
the two-dimensional hashing table EI(aq, a^) is sparse and further that each linked class only 
contains a few sites. Thus, the status-checking procedure should be very efficient. Indeed, 
our simulation implies that, close to p c , the CPU time for the status-checking procedure is 
of order 0(1) per site, irrespective of the dimensionality d. 

For site DP, one also sequentially visits the sites in queue Q t and generate via Eq. 
the next to-be-wetted sites one by one. However, the status- checking procedure is slightly 
different from that for bond DP. The rule is more restrictive. A newly calculated to-be- 
wetted site v t+ i cannot be added to queue Qt+i as long as it has already been visited, 

n 



irrespective of whether it was decided to be wetted or unwetted. Suppose v t+1 is derived 
by Eq. fl3]) from the kth element in queue Qt, one observes that, if vt+i has already been 
visited, it must be a neighbor of some earlier element in Q t with label k! < k. In other 
words, no element in Q t with label k' < k can be a neighboring site of v t+1 , if v 4+1 is to 
be wetted. This motivates us to classify the sites in Q t instead of Qt+i- Analogously, the 
classification is based on the first two coordinates of the elements in Q t , and all the elements 
in the same class are directly linked. When a to-be- wetted site v t+1 is calculated by Eq. fl3]), 
all the linked classes, of which the elements in Q t may be neighboring to v t +i, are identified. 
Within these linked classes, it is checked whether there exists an element in queue <Q) t which 
is neighboring to site v i+1 . If not, v t+1 is added on the top of queue Qt+i, and the hashing 
table and the linking are updated. By this way, albeit a potentially wet site v t+1 at t + 1 
may be visited several times and thus more than one random number is drawn, the status 
of v t+ i is determined by the first drawn random number. The status-checking procedure is 
only slightly more complicate than for the bond DP, and thus the efficiency remains. 

We mention that, near p c , the maximum computer memory needed in simulation is not 
determined by the average size of the wetted sites at time t max . Instead, the maximum 
length of the queues Q should be set as oc t^ ax , where the cut-off exponent y^- characterizes 
the wide-range probability distribution of the wet-site number M{t). Similarly, the hashing 
table M(xi,x 2 ) is of size oc ijj^ax, with the cut-off exponent y^ accounting for the probability 
distribution of the gyration radius 1Z of the wetted sites. It will be shown numerically that 
the cut-off exponents are in range 0.4733(<i = 1) < yN < l(d > 4) and 0.63267(<i = 1) > 
Ur > 1/2 (d > 4). In the actual implementation, it is monitored whether an overflow would 
occur in queues Q and the hashing table M. 

We conclude this section by mentioning that a different hashing technique, u sing 64-bit 



integers, was applied to study the standard percolation and the DP in Refs. 34j, |35 |. 



III. RESULT I: PERCOLATION THRESHOLD 

Using the algorithm in Sec. Ill Cj we study the bond and site DPs in (d+1) dimensions with 
d = 1-7. Preliminary and coarse simulations are performed for various p to approximately 
locate the percolation threshold p c by looking at the intersection of a dimensionless ratio Q^f 
for different system sizes, which will be defined later. Then, fine and extensive simulations 
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FIG. 4: Ratio Qjv versus p for some DPs. 
are carried out near the estimated critical point. 



A. Dimensionless ratio 

In the DP process, the number of wetted sites N(t) is recorded in computer memory 
for each time step t, and the statistical average N(t) = (N(t)) is calculated. An improved 
estimator of J\f(t) was used in Refs. 0, 35 1. and found to reduce the variance of N(t) by a 
factor of about 3 for some unknown reason. In this work, we did not apply this improved 
estimator. To reduce the statistical fluctuation, we simply measure the average m range 
(t, 2t] as 

(4) 



v 7 t'=t+i 



where t is taken as t = 2 l with i = 1, • • • , v ax — 1- Apparently, A/"(t) at the percolation 
threshold p c behaves the same as Af(t)-i.e., 77(t) oc t 9 . 

As discussed in the Introduction, in equilibrium statistical mechanics, universal dimen- 
sionless ratios are found to be very powerful in locating critical point. Therefore, we also 
define for DP a dimensionless ratio as 

Q N {t)=Jf{2t)/Jf{t). (5) 

According to the scaling behavior of A/"(£), we expect: (1), for p < p c , the Qj\r(p,t) value 
quickly approaches to as t increases, (2), for p > p c , it rapidly reaches 2 d , and (3), at 
p = Pc, it converges to the universal value 2 e . The Qj^f data versus p for different time 
steps t should display an intersection, which indicates the percolation threshold p c . This 
is indeed the case, and some results are shown in Fig. [4j A rough eye-view fitting yields 
p c = 0.644 70(1) for the square-lattice bond DP, p c = 0.382 222(1) for the d = 2 SC bond 
DP, p c = 0.207 918(1) for the d = 4 SC bond DP, and p c = 0.007 818 3(1) for the d = 7 
BCC bond DP. The estimate of p c for the d = 4 SC bond DP falsifies the existing result 
p c = 0.208 5(2). 
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TABLE I: Result from the separate fits of Qj\f in (1+1) dimensions. Superscripts b and s represent 
the bond and site DPs, respectively. 

According to renormalization group theory, the behavior of the wet-site number Af(p, t), 
and thus of the dimensionless ratio Q^f, obeys a universal scaling function near p c for all 



gb gs T b T s 

0.644 7001(2) 0.705 485 3(3) 0.478 025 0(4) 0.595 646 9(2) 



H b 



K b 



0.822 856 9(2) 0.839 9316(2) 0.658 968 9(2) 0.736 9317(2) 



TABLE II: Estimate of p c from the separate fits of Qn in (1+1) dimensions. 



DPs in the same universality class. In other words, we expect that 



Q*(p,t) = Q*{vt v *,ut v ») , 



(\P~Pc\ « Pc) 



(6) 



where v and u represent the amplitudes of the relevant and leading irrelevant scaling fields, 
respectively, and y p = 1/r/ii and y u < are the associated renormalization exponents. Func- 
tion Qj^f and exponents, y v and y u , are universal. The amplitude u characterizes the distance 
between the fixed point and the location of a critical system in the critical surface, and is 
thus model-dependent. The relevant scaling field v can be expanded in the occupation 
probability p, as 



where p C j is the percolation threshold for the jth DP model. The coefficients a, and ej are 
nonuniversal constants. 

In order to bring Eq. (ED in a suitable form to describe the Monte Carlo data for Qj^(p, t), 
the universal function can be Taylor-expanded as 



where the amplitude Qi, obtained from the ith derivative of Qn, is also universal. The 
universal value Q c relates to exponent 9 as Q c = 2 e . 

B. Separate analysis 

Result for d = 1. For d = 1, extensive simulations were carried out for the bond and 
site DPs on the square (S), triangular (T), honeycomb (H), and kagome (K) lattices, which 
include 8 models in total. The maximum time step is set as t max = 2 16 for all these 8 models, 
which corresponds to virtual lattices of ~ 10 9 sites. For each system, about 10 9 samples are 




(7) 



Qn = Qc + t Vp v J Q 1 + t 2y * Vj 2 Q 2 + ■■■ + u 3 t yu + ■■■ , 



(8) 



taken. The simulation runs on a cluster of 128 CPUs and takes about 12 years of CPU time 
for a PC of 2. 66 GHz frequency. 
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FIG. 5: Ratio Qj^ versus parameter x for DPs with d = 1. Parameter x is calculated as x = 
qij(p c j —Pj)t Vp , with qy and p c j taken from Tables. HI and ITT1 The excellent collapse indicates the 
universality class. 

According to a least-squared criterion, the Qn{p, t) data are fitted by 



which is obtained by substituting Eq. (J7J) in Eq. (jSJ) and ignoring high-order nonlinear and 
correction terms. For simplicity, the second correction exponent y<i is fixed at y^ = —2. In 
the fit, the Monte Carlo data for t < t m i n are excluded, where t m \ n is gradually increased 
until the residual x 2 is comparable to the degrees of freedom (d.f.) and the \ 2 value does not 
drops significantly when a few more data points are thrown away. Furthermore, analogous 
fitting procedures are carried out by excluding those poorly-determined terms or including 
additional terms in Eq. fl9]). Consistency and stability of the results are checked. The error 
margins are estimated by taking into account the statistical errors directly from the fits and 
the variation in different fits. The estimated error bars are taken so that they can at least 
cover the fitting results from different fits. 

The results are shown in Tables. HI and HT1 Universality is confirmed by the good agreement 
of the universal ratio Q c and the exponent y p for all the 8 models. We also compute 



Qm = Q c + QijiPcj - Pj)t yp + q2j(p C j - Pj) 

+ Cj ( Pcj - Pj )t V " +V " + feyt"" + b 2j t y2 , 



2 t 2Vp + q3 j {Pc j -p j ) 3 t 3y " 



(9) 



I I 



Model Q c 


V P 


Vu 




<72 


qi/q\ 


h 


sc£ 


1.1734(1) 


0.7762(5) 




-3.15(1) 2.65(5) 


0.267(7) 




sci 


1.1737(4) 0.7767(7) 


-0.46(8) 


-2.42(1) 1.56(1) 


0.266(4) 


-0.022(6) 


BCC^ 


1.1735(2) 0.7758(3) 




-4.20(1) 4.70(5) 


0.266(5) 




BCC S 2 


1.1729(4) 0.7761(6) 


-0.58(7) 


-2.92(2) 2.28(3) 


0.267(8) 0.027(5) 


SC^ 


1.0774(6) 0.905(4) 




-3.90(4) 6.4(2) 


0.42(3) 


— 


sci 


1.077(2) 


0.904(7) 


-0.33(5) 


-2.67(6) 3.2(2) 


0.45(5) 


0.028(5) 


BCC| 


1.0762(5) 0.902(5) 




-8.36(8) 30(1) 


0.43(3) 


— 


BCC1 


1.077(2) 


0.901(5) 


-0.37(3) 


-5.0(1) 


11.4(4) 


0.46(4) 


0.034(2) 


SC^ 


1.007(7) 


0.984(7) 


-0.20(6) 


-4.34(7) 10.4(6) 


0.54(5) 


0.031(6) 


SC| 


1.011(4) 


0.983(8) 


-0.28(7) 


-2.9(3) 


4.6(6) 


0.5(2) 


0.059(3) 


BCC^ 


1.007(3) 


0.984(3) 


-0.16(5) 


-15.0(3) 126(7) 


0.56(6) 


0.016(2) 


BCC| 


1.011(5) 


0.975(9) 


-0.30(6) 


-9.6(9) 


43(8) 


0.5(2) 


0.059(3) 


SC!? 


1.000(2) 


0.999(7) 


-0.45(9) 


-5.2(3) 


15(3) 


0.6(2) 


0.025(9) 


SC| 


1.000(2) 


0.997(5) 


-0.47(4) 


-4.2(1) 


11(2) 


0.6(2) 


0.052(5) 


BCC^ 


1.0000(3) 


1.001(3) 


-0.45(7) 


-30.9(7) 566(37) 


0.59(7) 


0.0075(7) 


BCC| 0.9994(6) 


1.005(9) 


-0.41(8) 


-21(2) 


235(33) 


0.5(2) 


0.04(1) 


sc£ 


0.9998(5) 


0.993(4) 


-0.8(3) 


-6.8(3) 


26(3) 


0.6(2) 


0.021(9) 


sci 


0.9999(5) 


0.998(4) 


-0.9(2) 


-5.4(2) 


16(1) 


0.55(8) 


0.05(1) 


BCC^ 


0.9999(4) 


0.999(5) 


-0.9(4) 


-64(2) 


2539(192) 


0.62(9) 


0.0018(5) 


BCQ 0.9999(5) 


1.003(4) 


-0.9(2) 


-48(2) 


1280(90) 


0.56(9) 


0.05(2) 


SC^ 


0.9997(3) 


1.000(4) 


-1.3(5) 


-7.6(2) 


33(3) 


0.57(9) 


0.014(5) 


SC S 7 


0.9996(7) 


1.003(4) 


-1.2(2) 


-6.4(2) 


24(2) 


0.59(9) 


0.040(8) 


BCC^ 


1.0000(2) 


1.002(4) 


-0.9(2) 


-126(4) 


10032(1344) 0.6(2) 


0.0005(2) 


BCC S 7 


1.0000(2) 


0.997(4) 


-1.1(1) 


-108(3) 


6698(900) 


0.6(1) 


0.034(8) 



TABLE III: Result from the separate fits of Qn for DPs with d = 2-7. Superscripts b and s are for 
the bond and site DPs, respectively, and subscript denotes the spatial dimensionality d. Symbol 
"— " means that the corresponding parameter cannot be well determined, since the error margin is 
comparable or larger than the absolute value of the parameter itself; thus, the corresponding term 
is not included in the fit. 

1 K 



d SC b 



BCC b BCC S 



2 0.382 224 66(6) 

3 0.268 356 38(10) 

4 0.207 918 20(7) 

5 0.170 615 20(4) 

6 0.145 089 92(4) 

7 0.126 387 50(4) 



0.435 314 3(1) 
0.303 395 52(9) 
0.231046 98(8) 
0.186 513 61(5) 
0.156 547 20(4) 
0.135 004 21(4) 



0.287 338 38(6) 
0.132 37419(3) 
0.063 763 404(7) 
0.031456 634(2) 
0.015 659 387(6) 
0.007 818 373(3) 



0.344 573 86(9) 
0.160 96131(10) 
0.075 58518(1) 
0.035 972 52(2) 
0.017 333 055(5) 
0.008 432 989(3) 



TABLE IV: Estimate of p c from the separate fits of Qj\f for DPs with d = 2-7. 

which is equal to Q2/Q\- It can be seen that the (fcAZi values are also identical for 
different systems, supporting the statement that Eq. ([6]) is universal. To further demonstrate 
the universality, we plot the Qjy{p,t) data for all these DPs versus parameter x, where 
x = QijiPcj — Pj)t Vp is calculated from Tables. [I] and [Til The excellent collapse is observed in 
Fig. |5l where the data for the small time step t < 1536 are excluded in order to avoid the 
influence from finite-size corrections. 

The fits for all the 8 models give the leading correction exponent as y u ~ —1. The 
estimates of p c are already quite precise, of which the error margins are at the 7th decimal 
place. 

Result for d = 2 and 3. The simulation for d — 1 and 2 did not use the linking-hashing 
technique, or the cumulative-probability trick. For d > 2 these two tricks improves the 
efficiency and limits the used computer memory. For the BCC lattice in (7+ 1) dimensions, 
the coordination number is A = 2 d = 128; no significant decreases of the efficiency is 
observed. 

There exist 24 systems for the site and bond DPs on the SC and BCC lattices for d = 2- 
7. For the majority of these systems, the maximum time step is set at t max = 2 13 ; for the 
others, t max is larger than 2 13 . For each of these systems, about 10 9 samples are taken. The 
total CPU time is about 43 years for a PC with 2.66GHz frequency. 

Following the similar procedure for d — 1, the data are fitted by Eq. 0. The 
results are shown in Tables. inil and lTVl The determinations of Q c , y p , and q2j/q\j agree well 
with each other for all the DPs in the same spatial dimension. Moreover, in two and three 
dimensions, the Q c = 2 e value clearly deviates from 1, implying the critical exponent 8^0. 

1 a 



The estimate of p c is rather good, the precision of which already exceeds that of the existing 
results. However, it is rather difficult to determine the leading correction exponent y u . For 
the bond DP, the Qj^ data can be described by Eq. (Q without including the correction 
term with exponent y u \ this is represented by symbol "-" in Table IIHI If such a term is 
included, the amplitude b\ and the associated error margin are of the same magnitude. For 
the site DP, we determine the correction exponent as y u ~ — 0.50(g? = 2) and — 0.35(d = 3). 

Result for d = 4. The upper critical dimensionality for the DP universality class is 
believed to be d c = 4. For d > d c , the mean-field prediction becomes valid and the critical 
behavior is governed by the Gaussian fixed point in renormalization group theory, which 
yields the renormalization exponent 1/i/m = 1 and 6 = 0. Nevertheless, at d = d c = 4, 
the existence of dangerous irrelevant scaling fields typically leads to both multiplicative and 
additive logarithmic corrections for finite-size scaling. 

Despite of the possible existence of logarithmic correction, as a crude analysis we also 
fitted the Q_\f data by Eq. (Q. The results are shown in Tables [TTT1 and HVl While slowly 
convergent corrections with exponent y u « —0.2 are observed in most cases, the data 
for the bond DP on the BCC lattice do not suffer from severe finite-size corrections. The 
estimates of Q c and y p agree more or less with the mean-field prediction Q c = 1 and y p = 
= 1, albeit significant deviations can indeed observed for some cases. From these 
deviations, it seems that multiplicative logarithmic corrections with negative exponent exists 
in terms associated with y p . 

For d = d c , the leading and subleading, logarithmic corrections have been calculated by 
Janssen and Stenull from field theory |36|, and are partially confirmed by Monte Carlo 



simulation 



35]. At critical point, the scaling behavior of the average wet-site number J\f(t) 



is described by |36| 

Af(t) = N [ln — - Mnln — + a]° (10) 

to t\ 

where the known quantities are a = 1/6, b = 1.30204, and a = 0.1831, and unknown 



constants No, t and t\ are model-dependent. It is further predicted [37| that, if only the 



leading logarithmic corrections are taken into account, J\f(p,t) behaves as 

N(p,t) w (\nt) ai $x(t(\nt) a2 Ap) (11) 
where $jv" is a universal scaling function and the exponents are cti = 1/6 and «2 = —1/6. 



i - 



Accordingly, we fit the J\f(p, t) data by 
Af{p,t) ={\nt + h 1 ) ai x 



(12) 



k=l 

with yi — — 1 and 2/2 = —2 fixed. Indeed, the A/"(p, t) data can be well described by Eq. ffl~2]) 
with «i = 1/6 and «2 = —1/6 being fixed; the results for hi and /12 are shown in Table IVl 
If hi and h 2 are fixed at the values in Table IVl the exponents ai and a 2 can be determined 



Model SC^ SC| BCC^ BCC| 
hi 0.6(2) -1.27(20) 3.4(2) -1.1(2) 
h 2 2.1(6) -2.2(7) 3.2(6) -2.0(5) 



TABLE V: Estimate of hi and hi in Eq. (fl~2j) . with ai = 1/6 and 02 = —1/6 being fixed. 

as in Table IVIl in agreement with the theoretical values. Nevertheless, the present Monte 
Carlo data do not allow the simultaneous determination of ct^ and hi (i = 1,2). 



Model SC h 4 SC| BCC^ BCQ Ref [37] 
oji 0.171(8) 0.167(4) 0.165(3) 0.166(5) 1/6 
a 2 -0.17(2) -0.17(2) -0.16(2) -0.16(2) -1/6 

TABLE VI: Estimate of ai and a 2 in Eq. (]12p . with hi and h 2 fixed at the values in Table IVl 

We also fit the Qj\f(p, t) data by 

/ 1 9 \ ai 

Qn = Qc ( 1 + ; " r + c,-(p c -p)^ + ^(lnt + /i 2 ) a2 



In t + /i! 

3 

+ <?fc(Pe - p)¥^(lnt + /i 2 ) fcQ2 + M V1 + M y2 , (13) 
fc=i 

where ai = 1/6, 02 = — 1/6, ?/i = — 1, y 2 = — 2 are fixed. To estimate Q c and y p , we took the 
values of hi and h 2 in Table. |V] The results are shown in Table [VTTI the estimates of Q c and y p 
are now consistent with the mean-field predictions. With Q c = 1 and y p = 1 being fixed, the 
estimate of the percolation threshold is p c (SC^) = 0.207918 18(4), p c (SC 8 4 ) = 0.231046 87(4), 
p c {BCC\) = 0.063 763 395(4), andp c (BCQ) = 0.075 585 16(2), in good agreement with those 
in Table El 



Model Q c y p 


Ql Q2 


Q2/qj 


bi 


SC^ 


1.0002(4) 1.002(7) 


-5.50(8) 16(2) 


0.53(9) 


-0.15(4) 


SQ 


1.0000(5) 1.004(4) 


-3.48(6) 6.4(4) 


0.52(6) 


-0.6(2) 


BCC^ 


0.9999(2) 1.000(2) 


-19.6(2) 216(10) 0.56(4) 


-0.06(4) 


BCC| 


1.0003(7) 1.003(3) 


-10.3(2) 55(4) 


0.52(6) 


-0.05(2) 



TABLE VII: Result of the separate analysis for (1 = 4, according to Eq. (|13|) . 

Result for d = 5, 6, and 7. For d > d c = 4, the critical behavior of DP should be 
described by mean-field theory, which yields the renormalization exponent y v — 1/i/n = 1 
and 8 = 0. This implies that, on average, the wet-site number J\f(t) oc t e = 1 is just of a few 
sites, independent of the propagation time. It will be shown later that, at d — d c , one has 
Af(t) « 1.2 up to t = 10 4 . 

The data are fitted by Eq. (jHj), and the results are shown in Tables II III and IIV1 
The mean-field predictions Q c = 2 e = 1 and y p = 1 are well confirmed. The estimated 
values q\ijq\ are identical for all the site and bond DPs on different lattices and in different 
dimensions. The leading correction exponent is y u ~ — 1/2 for d = 5 and y u ~ — 1 for d = 6 
and 7. We suspect that = 5) ~ — 1/2 is due to the crossover phenomena arising from the 
small distance to the upper critical dimensionality. The determination of the percolation 
threshold is up to the 8th or even the 9th decimal place, reflecting the efficiency of our 
algorithm. 

C. Simultaneous analysis 

The results of Sec. IIIIBI indicate that, as expected, the FSS behavior of observables near 
the critical point is indeed described by universal scaling function for systems in the same 
universality class. For the present work, this means that all the DPs in the same dimension 
d are in the same universality class; the number of different systems is 8 for d = 1 and 4 for 
d > 1. Assuming that the universality strictly holds and thus the universal parameters are 
exactly the same, we combine the Qn{p, t) data for all the DPs with the same d and analyze 
them simultaneously by a single formula, in which the universal parameters only occur once. 

i n 



Separating the universal from the nonuniversal constants in Eqs. (J7J) and (jSJ) yields 

Qn = Qc + Qiajipcj - Pj)t yp + Q 2 a 2 j ( Pcj - Pj)H 2v * + Qstfipcj - Pj )H 3Vp 

+ cjipcj - Vj )t^ + b Xj t*> + b 2j t y * , (14) 
where Q c , Qi (i — 1, 2, 3), y p , and y u are universal. 





d Q c Q2 Q3 






1 1.24293(5) -0.33(1) -0.29(1) 


0.5766(2) -1.03(3) 




2 1.17346(13) 0.266(2) -0.041(8) 0.7758(4) -0.52(4) 




3 1.0769(2) 0.441(4) 0.09(3) 


0.9033(16) -0.39(2) 




4 1.0075(7) 0.544(8) 0.24(3) 


0.983(3) -0.22(3) 




5 1 (fixed) 0.56(2) 


1 (fixed) -0.48(2) 




6 1 (fixed) 0.56(2) 


1 (fixed) -0.95(8) 




7 1 (fixed) 0.57(2) 0.37(9) 


1 (fixed) -1.12(6) 


TABLE VIII: Result for the universal parameters from the simultaneous fit of Qj^. 


Lattice 


Site 


Bond 




p c (Present) p c (Known) 


p c (Present) p c (Known) 


square 


0.705 485 25(9) 0.705 485 22(4) [20J 


0.644 70019(5) 0.644 700185(5) [20J 




0.705 489(4) [21] 


0.644 70015(5) [19] 


triangular 


0.595 646 86(8) 0.595 646 75(10) [22] 


0.478 025 26(9) 0.478 025 25(5) [22] 




0.595 646 8(5) [19] 


0.478 025(1) [12] 


honeycomb 


0.839 93170(5) 0.839 933(5) [18] 


0.822 856 90(4) 0.822 856 80(6) [22] 


kagome 


0.736 93175(8) 0.736 93182(4) [22J 


0.658 969 03(7) 0.658 96910(8) [22] 



TABLE IX: Estimated percolation threshold p c of the DPs in (1+1) dimensions, obtained from the 
simultaneous fit of Qj\f- 

Following the similar procedure as in Sec. IIIIB[ the Qj\r(p,t) data near p c are fitted by 
Eq. (JT4"]) according to the least-squared criterion. Note that, in Eq. (fl4|) . the universal 
parameter Qi and the model-dependent constants a,j are not independent parameters and 

on 



Lattice 


Site 

p c (Present) p c (Known) 


Bond 

p c (Present) p c (Known) 


(2+l)-d SC 


0.435 31419(6) 0.435 31(7) [24] 


0.382 224 64(4) 0.382 223(7) [24] 


(2+l)-d BCC 


0.344 573 98(6) 0.344 573 6(3) [28] 
0.344 575(15) [27] 


0.287 338 38(4) 0.287 338 3(1) [25] 
0.287 338(3) [24] 


(3+l)-d SC 


0.303 395 54(5) 0.302 5(10) [26] 


0.268 356 32(4) 0.268 2(2) [14] 


(3+l)-d BCC 


0.160 96123(3) 0.160 950(30) [27] 


0.132 374 21(2) not found 


(4+l)-d SC 


0.231 046 94(3) not found 


0.207 91819(2) 0.208 5(2) [14] 


(4+l)-d BCC 


0.075 585173(6) 0.075 585 0(3) [28] 
0.075 582(17) [27] 


0.063 763 402(3) not found 


(5+l)-d SC 


0.186 513 57(4) not found 


0.170 61518(3) 0.1714(1) [14] 


(5+l)-d BCC 


0.035 972 547(8) 0.035 967(23) [27] 


0.031456 634(1) not found 


(6+l)-d SC 


0.156 547 22(2) not found 


0.145 089 95(3) 0.145 8 [14] 


(6+l)-d BCC 


0.017 333 055(3) not found 


0.015 659 385(3) not found 


(7+l)-d SC 


0.135 004 21(2) not found 


0.12638751(2) 0.1270(1) [14] 


(7+l)-d BCC 


0.008 432 989(2) not found 


0.007 818 373 0(8) not found 



TABLE X: Estimated percolation threshold p c of the DPs in (d+1) dimensions with d = 2-7, 
obtained from the simultaneous analysis of Qj\f. 

thus cannot be both determined in the analysis. We set Qi = 1 in the fit, which defines a 
scale of the relevant field. The results are summarized in Tables IVIIIllXl For DPs above the 
upper critical dimensionality d c = 4, the mean-field predictions, Q c = 1 and y p = 1, are well 
confirmed if Q c and y p are left as free parameters to be determined by the fit. Thus, we only 
present in Tables IVIIIllXl the results with Q c — 1 and y p = 1 being fixed. It is suggested 
that the Q2 value is identical to a single value for d > d c , and thus the data for all 
the DPs for d > d c can be in principle combined together and analyzed simultaneously. For 
simplicity, we did not carry out such an analysis. 

At d c = 4, the results Q c = 1.0075(7) and y p = 0.983(3) clearly deviate from the mean- 
field values Q c = 1 and y p = 1, due to logarithmic corrections. Thus, we fit the data 



01 



by 

3 

+ Y,^a){p cj -p j ) k t ky -{\nt + h 2j T" + b lj t^+b 2j t^ , (15) 

k=l 

where a± = 1/6, a 2 = —1/6, y\ = —1, y 2 = —2, Q c = 1 and y, p = 1 are fixed. The 
determined percolation thresholds, p c (SC\) = 0.20791818(3), p c (SC£) = 0.231046 88(4), 
p c (BCC^) = 0.063 763 395(3), p c (BCC|) = 0.075 585 15(1), agree well with those in Table D3 
Compared to the separate analysis, the precision of p c and of universal parameters is 
significantly improved. A comparison of the threshold p c with the existing results is presented 
in Tables ILXl and 1X1 In (1+1) dimension, our determination of p c agrees with the existing 
results, with a comparable or somewhat better precision. These d = 1 results are mostly 
obtained by series expansions. We consider it valuable to provide an independent study by 
a distinct method. For d > 2, Monte Carlo simulation turns out to be an efficient method. 
Our results for d = 2-7 provide a number of first estimates of p c or a significant improvement 
on the precision of p c . 

Finally, let us discuss the p c values given in Table 1X1 for d > 2. In terms of coordination 
number A, Kurrer et al. suggested an empirical formula 



l/p c ^ai + a 2 A, (16) 

where (i=l, 2) are some unknown constants. The formula (fT5"|) applies for large A, and 
cannot describe the p c values for small A. In Fig. [6j we plot the l/p c value versus A. It is 
interesting to observe that the slopes are approximately identical for the bond and site DPs 
on the SC lattice, while on the BCC lattice they clearly have different slopes. We take the 
d = 6 and 7 data, and calculate a\ and a 2 according to Eq. (fi~6l) . The results are shown in 
Table IXll On this basis, we suspect: (1), the a 2 value is not necessary equal to 1, and (2), 



sc b 


SC S 


BCC b 


BCC S 


d -0.2470 
a 2 1.0199 


-0.7474 
1.0193 


-0.1849 
1.0007 


-3.1954 
0.9514 



TABLE XI: Estimate of a\ and a 2 in Eq. f)16[) . calculated from the d = 6 and 7 data, 
the a 2 value is identical for the bond and the site DP on the SC lattice. 




FIG. 6: 1/pc versus coordination number A for the bond and the site DP on the BCC lattice. The 
lines are obtained by Eq. (|16p according to the d = 6 and 7 data. The inset displays for the SC 
lattice. 

IV. RESULT II: SCALING AT CRITIC ALITY 



Sampled quantities and scaling behavior. To study the scaling behavior of other quantities, 
we simulate the bond DP on square and higher-dimensional BCC lattices with d from 2 to 5, 
at p c = 0.644 700 185(d = 1), 0.287338 38(d = 2), 0.132 374 21(^ = 3), 0.063 763 404(^ = 4), 
and 0.031 456 633(<i = 5), which are, within the error margin, at the estimated critical points 
in Tables HXl and IXl Besides the wet-site number N(t), we sample the wetting probability 
P(t) and a revised gyration radius R(t). The wetting probability P(t) accounts for the 
probability that, at time t, the DP propagation still survives. Namely, P(t) = 1 represents 
that there is at least one wetting site at time t, while P(t) = is for no wetting site. The 
revised gyration radius R(t) is defined as 



R(t) 







N(t) = 



with R s (tf 



N(t) 

Eiv. 

i=i 



(17) 



y/R s (t) 2 /N(t) N(t) > 1 
where |vj| is the modulus of the spatial coordinates of site i. On the basis of P (t), N(t) 



no 



and R(t), we measure 

V(t) = (P(t)) 

Af(t) = (N(t)) , A4(t) = (N(t) k ) , 

K(t) = (R(t)) , K k (t) = (R(t) k ) , (18) 
with k = 2,3,4. Note that TZ(t) is different from the traditional definition 

n T {tf = (R s (tf)/(N(t)) = n 2 (t)/v(t) . (19) 

In Monte Carlo simulation, the revised gyration radius TZ{t) is a directly measured observ- 
able, while the traditional radius TZrit) is a composite quantity. 

At critical point p c , it is already known that the scaling behavior of V(t), Af(t), and 
K T (t) is 

V(t) oc t~ 5 , M(t) oc t 9 , and K T {t) oc t 1/z , (20) 

where critical exponent z characterizes the anisotropy between the temporal and the spatial 
direction. Let z/ii and u± denote the exponents describing the behavior of the correlation 
lengths along the temporal and the spatial directions, respectively, £y oc \p — Pc| _!/|1 and 
oc \p — p c \~ u± , one has z = v^jv^. From Eqs. ffTTj) and f|T9|) . it can be derived that, at 
criticality, the revised gyration radius TZ{t) behaves as 

K(t) oc V(t)K T (t) oc r s+1/z . (21) 

Analogously, it is expected that, at criticality, the higher-order moments of the wet-site 
number N(t) and of the revised gyration radius R(t) behave as 

J\f k (t) oc t m k and K k (t) oc t v *k , (k = 2,3,4). (22) 

Probability distribution. We further examine the probability distribution of N(t) and 
R(t). For simplicity, we consider the conditional probability that there is at least one wet 
site. More precisely, we define F^(N = s,t) as the probability that, conditioning on that 
the DP propagation still survives, N(t) is equal to given number s. The same applies to the 
conditional probability Pr(-R = s,t) of R(t). 

According to FSS theory, together with the Monte Carlo data, we assume that the prob- 
ability distributions, P(iV = s,t) and F(R = s,t), are described by 

F N (N = s,t)<x r yN F N (s/t yN ) , F R (R = s,t)<x r y »F R ( s /t y *) , (23) 




s/t yN s/t yR 



(c)Pjv for d = 5 (f)F R for d = 5 

FIG. 7: Probability distributions Pat and P#. The exponents yjv an d 2/_r are taken from Table IXlIl 
The times steps are t = 2048, 4096, •• • ,131072 for d=l, t = 512, 1024, •• • ,32768 for d=4, and 
t = 256,512,- •• ,16384 for d=5. The slope of the dashed straight line is 1/un — 1 for P^ and 
1/Vr - 1 + d for F R . 
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each involves only one critical exponent, y^ or y R . Functions, F N and F R , are expected 
to be universal. Figure H displays F N t VN {¥ R t VR ) versus s/t VN {s/t VR ) for d = 1, 4 and 
5. The exponents are set to be y^ = 0.4733 and y R = 0.63267 for d — 1, as numerically 
determined latter. For d = 4 and 5, the mean-field predictions y^ — 1 and y R = 1/2 are 
used. A good collapse is observed. Theoretically, for d = d c = 4, Eq. (1211 suggests a better 
collapse should occur for F N tY-s ~ s/(tY~s) and P^y 24 ~ s I (t^Y 24) . Nevertheless, no 
significant improvement is observed in the eye view. 

The dashed lines in Fig. [7] have slope — 1 for P^ and l/y R — 1 + d for F R . The 

approximate collapse of the line and the numerical data indicates that for small s, the 
universal functions and F R behave as a power law, as typical in critical phenomena. 
Accordingly, it seems that one can rewrite Eq. f[2"3"|) as 

¥ N (N = s,t) oc r l s llvN - l f N (s/t VN ) 

F R (R = s,t) oc t- l - dy *s l ly*- l+d f R (s/t y *) , (24) 

where fN(s/t VN ) and fti(s/t yR ) are universal. 

From the probability distribution (123]) . it can be derived that 

M{t) oc r s+VN , M k {t) cxt~ 5+kyN 

K(t) oc r s+VR , ^ fc (t) oc r 5+kVR . (25) 

Comparing Eq. ( 1251) to Eqs. ( 1201) and ( 12TT) . one obtains the relations 

y N = 9 + 5 , and y R = l/z . (26) 

As well known, for DP there exist only three independent exponents, which are usually 
chosen as 0, v\\, and u±. The others can be obtained by scaling relations js] 

S = /3/ U[h z = u ^ u±i and e = (du ± -2p)/v\\, (27) 

where the last one involves the spatial dimensionality d and is usually referred to as the 
hyper-scaling relation. In this work, for the convenience of finite-size scaling, we choose y p , 
y R and S as the three independent exponents. Then, the other exponents can be obtained 

as 

^|| = l/y P , z = l/y R , v L = y R /y p , 

9 = dy R - 26 , (5 = 5/y p , y N = dy R -5. (28) 



At d = d c — 4, the scaling behavior (|25|) is modified by multiplicative logarithmic fac 



361 ] . it can be derived that a power-law 



tor. From the calculations by Janssen and Stenull 
behavior normally becomes 

t x i*-- f -pn— - 6 In In - + a,] 01 - := f kf T ai , (29) 
to tl 

where x is some critical exponent, and x m .f. represents the mean- field value. Amplitude b is 
universal and equal to b = 1.302 04, i an d £i are unknown integration constants from the 
renormalization-group flow, which are model-dependent but observable- independent, and a x 
is observable-dependent. One further knows that 

5 m s. = 1, a s = -1/2 , and m . f . = , ate = 1/6 , 

J/jv.m-f. = 1 , « yjv = - 1 / 3 , and yR,m.i. = 1/2 , a te = 1/24 . (30) 

On this basis, we have 

V{t) = V t- l Y 1 ' 2 , 

M(t) = Af Q t°Y 1/6 , N k (t) = tf kQ t k - 1 Y<- 3 - 2k V 6 , 

Hit) = TZ Q t- l ' 2 Y lz ' 2i , TZ k {t) = n k0 t^ 2 Y^ 12+k ^ . (31) 

In order to numerically confirm Eq. ( 13T1) . we plot the critical J\f(t), V ■ t, and 712 data for 
the d = 4 BCC bond DP in Fig. ED^a-c), respectively; the leading behavior of these quantities 



is expected by the 
-1.5193 for V-t 



ogarithmic factor. Amplitude a in Eq. (129]) is known as 0.1831 for A/", 
]. With a = 0.1831 and ag — 1/6 being fixed, the analysis of the Af data 
yields to ~ 0.0072, t\ ~ 0.037 and iVo ~ 0.8119. Since to and t\ are observable-independent 
for given system, they are fixed in the analysis of the V ■ t and IZ2 data. The results are 
shown in Fig. [S(a-c). For comparison, the data at p — 0.063 763 318 and 0.063 763 490 are 
also shown. The clear deviation from the theoretical prediction for large t implies that 
they are indeed respectively in the dilute and the absorbing phase in the DP process. One 
also expects that, from Eq. (l3T]) . the logarithmic factors Y in product NNz/t would cancel 
out and thus AOV^/t rapidly approaches to a constant as t increases. In Fig. [SJd), we 
plot A/A/2/i — O.llt -0,25 versus t, where — O.llt" 025 is obtained from the fit and the small 
exponent —0.25 is probably due to logarithmic correction. 

Determination of critical exponents. From the Monte Carlo data at criticality, we aim to 
numerically determine the critical exponents S, 9, z, y^, yn, etc. Besides V(t), A/"(i), and 
K(t), we also analyze Af k (t), K k (t), Af k (t)/M k -i(t) and ll k {t) /TZ k ^{t) with k>2. 
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FIG. 8: 

in (a-c) 
dashed 



Plot of various quantities for the (4+l)-dimensional BCC bond DP. The analytical curves 
are calculated from Eq. (|3ip with the unknown parameter obtained from the fits. The 
line in (d) is just for guiding eyes. 



To reduce one more unknown parameter in the fitting formula, we measure dimensionless 
ratios as following 



Q P (t) 

QN k 

D Nk 

which are fitted by 



Q 



Rk 



V(t) 
N k (2t) 

A4(t) 
A4(2t)A4-i(t) 
A4(t)A4_i(2t) 



M{t) 
1l k (2t) 

n k (t) 



nit) 



D 



Rk 



n k {2t)n k ^{t) 

TZ k (t)TZ k ^(2t) 



(fc = 2,3,4) 



(32) 



Q = Q C + ht^ + b 2 t^ + b 3 t y3 , (33) 
with y 2 = —2, y 3 = —3. In comparison with Eq. (Q, all the p-dependent terms are absent 



since it is at critical point. For d = d c , taking into account the logarithmic corrections, we 
have the formula as 

Q = &(! + rr^rT + h ^ vi + h * tm + b ^ tm > ( 34 ) 

mt + n 

where a can be obtained from Eq. (l3Tj) and y\ = —1. 

From the estimate of Q c , we calculate the associated exponent from the relation Q c = 2 X . 
For example, the fit of the Q-p(d = 1) data gives Q c ,v = 0.895 37(1), which yields 5 = 
— In Q Cy p/ In 2 = 0.159 44(2). The same applies to other quantities. The results are shown 
in Tables IXHl and IXHIt the latter is for the 5 result. The symbol "*" for d = 4 represents 
that the results are from the fits with logarithmic correction. 

All the exponents shown in Table 1X111 are mutually self-consistent. Further, the scaling 
relation, Eq. ( )26l) . is well confirmed. For d < d c = 4, the hyper-scaling relation in Eq. ( 1281) 
also holds. 

On the basis of Tables. IVIIII and IXII| we calculate other critical exponents from the 
scaling relations ( 128]) . Note that we directly take y^ — 8 in Table IXIII as our estimate of 
exponent 9 instead of calculating it by 9 = dy^ — 25, avoiding further error propagation. 
The estimates of these critical exponents, together with the existing results, are shown in 
Table IXIIII For d — 1 , our results are consistent with the existing result but worse than the 
latter. For d = 2, 3, our estimate significantly improves over the existing results by a factor 
from 4 to 10. For d — d c — 4, the results with the logarithmic corrections agree well with 
the mean-field predictions. For d = 5, the mean-field values are confirmed. 

V. CONCLUSIONS 

We present a systematic and high-precision Monte Carlo study of the bond and site di- 
rected percolations on (d + l)-dimensional lattices with d from 1 to 7. In order to obtain 
accurate results, the simulation of large time-steps is necessary. However, due to the restric- 
tion of computer memory, the large-time-step simulation is notoriously difficult. We formu- 
late an efficient algorithm by adopting the cumulative-probability and the linking- hashing 
technique. This allows us to obtain high-precision Monte-Carlo data within reasonable com- 
puting source. 

A dimensionless ratio Qj\f is defined and used to determine the percolation threshold 

on 



d+1 




2yjv-<5 3y N -5 


Vn(Dn 2 ) 


Vn{D Nz ) 


yN(D Ni ) 




Vr-8 


^Vr-S 3y R -5 Ay R -5 


Vr{Dr 2 ) 


yR(D R:i ) 


yR.(D R4 ) 


1+1 


0.31365(7) 0.7870(2) 1.2603(4) 1.7338(6) 0.4733(2) 


0.4734(2) 


0.4735(3) 




0.47313(4) 1.1057(2) 1.7386(3) 2.3714(2) 0.63267(6) 0.63275(8) 0.63274(8) 


2+1 


0.23081(7) 0.9123(5) 1.5937(9) 2.275(2) 


0.6815(4) 


0.6814(5) 


0.6816(7) 




0.1154(6) 


0.6813(3) 1.2472(3) 1.8133(3) 0.56607(6) 0.56608(7) 0.56610(6) 


3+1 


0.1073(8) 


0.9524(6) 1.7974(5) 2.643(1) 


0.8449(6) 


0.8445(9) 


0.845(1) 




-0.2102(5) 0.3160(7) 0.8426(5) 1.369(2) 


0.5271(2) 


0.5270(3) 


0.5269(4) 


4+1 


0.012(3) 


0.992(6) 1.962(3) 2.95(2) 


0.980(7) 


0.978(6) 


0.976(6) 




-0.456(4) 


0.049(4) 0.553(4) 1.059(7) 


0.5039(4) 


0.5039(5) 


0.5036(8) 


4+1* 


0.002(2) 


1.000(3) 2.01(2) 3.02(3) 


1.005(8) 


1.01(1) 


1.003(7) 




-0.494(3) 


0.006(5) 0.506(4) 1.007(5) 


0.5005(3) 


0.5000(3) 


0.4996(4) 


5+1 


0.0003(6) 


0.9994(6) 2.0002(4) 2.999(3) 


0.9992(6) 


0.9991(8) 


1.000(2) 




-0.498(1) 


0.002(2) 0.502(2) 1.002(2) 


0.4996(6) 


0.4996(7) 


0.4996(6) 



TABLE XII: Result for critical exponents from the dimensionless ratios in Eq. (|32p . Hereby, 
VA{D Ak ) (A=N, R; k=2, 3, 4) is from the fit of D Ak . The symbol "*" means that the loga- 
rithmic corrections are taken into account for d = 4. 

p c . In comparison with the existing results, our estimates about p c is comparable or some- 
what better for d = 1 and is much more precise for d > 1 dimensions. For many models, 
the p c estimates are presented for the first time. This suggests that, as equilibrium case, 
dimensionless ratios are also useful in locating the transitions in non-equilibrium systems. 

Based on the numerical data, the conditional probability distribution of the wet-site 
number N(t) and of the revised gyration radius R(t) at criticality is conjectured, which 
only involves one critical exponent. By means of finite-size scaling, we determine the critical 
exponents to a high precision, which can serve as a testing ground for future study of systems 
in the directed-percolation universality class. 



on 



d+1 


P 


v \\ 




Z 


9 


5 


l+l(Present) 


0.2765(2) 


1.7343(6) 


1.0972(6) 


1.5806(2) 


0.31365(7) 


0.15944(2) 


Ref. [20] 


0.276486(8) 1.733847(6) 1.096854(4) 1.580745(10) 0.313686(8) 0.159464(6) 


2+1 (Present) 


0.5812(6) 


1.2890(7) 


0.7297(5) 


1.7665(2) 


0.23081(7) 


0.4509(2) 


Ref. [24, 39J 


0.584(4) 


1.295(6) 


0.734(4) 


1.76(3) 


0.230 


0.451 


Ref. [25] 








1.7666(10) 


0.2303(4) 


0.4509(5) 


3+1 (Present) 


0.817(2) 


1.107(2) 


0.583(2) 


1.8972(8) 


0.1073(8) 


0.7382(5) 


Ref. [40J 


0.81(1) 


1.105(5) 


0.581(5) 


1.90(1) 


0.12 


0.73 


4+1 (Present) 


0.976(8) 


1.017(3) 


0.512(2) 


1.985(2) 


0.012(3) 


0.960(5) 


4+l(Present)* 0.996(3) 


1 (fixed) 


0.5005(5) 


1.998(2) 


0.002(2) 


0.996(3) 


5+1 (Present) 


0.9980(9) 


1 (fixed) 


0.4998(8) 


2.001(3) 


0.0003(6) 


0.9980(9) 


M.F. 


1 


1 


i 

2 


2 





1 



TABLE XIII: Final estimate of critical exponents in 1+1 to 5+1 dimensions. 
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